clc;clear;clearvars;
addpath('CEC2008\');
global initial_flag

% num_vari维度,max_evalutions最大评估次数,num_initial初始化点的个数,evalutions已经评估个数
num_vari = 100;
%max_evalutions = 5000 * num_vari;
max_evalutions = 500 * num_vari;
num_initial = num_vari;

% 搜索区间:f1[-100,100],f2[-100, 100],f3[-100,100],f4[-5,5],f5[-600,600],f6[-32,32],f7[-1,1],每一维搜索区间一样
upper_bound = [100,100,100,5,600,32,1];
lower_bound = [-100,-100,-100,-5,-600,-32,-1];

for func_num = 1:7

    % 使initial_flag==0满足此条件,换一个函数initial_flag重置为0
    initial_flag = 0;
    evalutions = 0;
    varphi = 0; %100维0,500维0.1和0.05,1000维0.15和0.1

    % 生成num_initial个种群,并获得其评估值
    sample_x = lhsdesign(num_initial, num_vari).*(upper_bound(func_num) - lower_bound(func_num)) + lower_bound(func_num).*ones(num_initial, num_vari);
    sample_y = benchmark_func(sample_x,func_num);
    evalutions = evalutions + size(sample_y,1);
    v = zeros(num_initial,num_vari);  %v初始化为0
    fmin = min(sample_y);
    arr_fmin = [];
    arr_fmin = [arr_fmin;fmin];

    while evalutions <= max_evalutions

        % 根据种群熵来划分pb和pw,种群熵的计算与MATLAB中entropy计算熵类似
        n_bin = num_initial;    %n_bin表示区间数
        res = rescale(sample_y);    % rescale将数组的条目缩放到区间 [0,1]
        p = imhist(res,n_bin);  % rescale将数组的条目缩放到区间 [0,1]
        p(p==0) = [];   % 除去p中的0
        p = p ./ numel(sample_y);% 正则化p使得sum(p)为1,numel计算数组中元素的数目

        % 为了E最大为1,所以用logn(p),而不是log2(p)
        n_nozero = size(p,1);
        E = -sum(p.*(log(p)./log(n_nozero)));
        
        % 划分pb和pw
        [sample_y,ind]= sort(sample_y); %按升序对 A 的元素进行排序
        sample_x = sample_x(ind,:);
        v = v(ind,:);
        mean_x = mean(sample_x);
        d = 0.25;   %100维0.25,500维0.35,1000维0.45
        if E > 1 - d
            mb = -(num_initial / d) * E + num_initial / d;
        else
            mb = num_initial;
        end
        mb = ceil(mb);
        mw = num_initial - mb;
        pb = sample_x(1:mb,:);  %越小越好
        pby = sample_y(1:mb,:);
        pbv = v(1:mb,:);
        pw = sample_x((mb + 1):num_initial,:);
        pwv = v((mb + 1):num_initial,:);

        next_samplex = [];
        next_sampley = [];
        next_v = [];
        % 更新pw
        for i = 1:size(pw,1)
            if size(pb,1) == 0  %pb有为空的问题
                xj = sample_x(1,:);   %当前最小值
            else
                xj = pb(randi(size(pb,1)),:);
            end
            pwv(i,:) = rand(1, num_vari) .* pwv(i,:) + rand(1, num_vari) .* (xj - pw(i,:)) + varphi .* rand(1, num_vari) .* (mean_x - pw(i,:));
            xi = pw(i,:) + pwv(i,:);
            xi(xi > upper_bound(func_num)) = upper_bound(func_num);
            xi(xi < lower_bound(func_num)) = lower_bound(func_num);% 范围检查
            if evalutions <= max_evalutions
                yi = benchmark_func(xi,func_num);
                if yi < fmin
                    fmin = yi;
                end
                arr_fmin = [arr_fmin;fmin];
                fprintf("func_num:%d evalutions:%d fmin: %.4f\n", func_num, evalutions, fmin);
                evalutions = evalutions + 1;
                next_samplex = [next_samplex;xi];   %加到下一次迭代中
                next_sampley = [next_sampley;yi];
                next_v = [next_v;pwv(i,:)];
            else
                break;
            end
        end

        % 更新pb
        while ~isempty(pb)
            if size(pb,1) == 1
                next_samplex = [next_samplex;pb];
                next_sampley = [next_sampley;pby];
                next_v = [next_v;pbv];
                pb = [];
                pby = [];
                pbv = [];
            else
                k1 = randi(size(pb, 1));    %随机抽取两个点
                k2 = randi(size(pb, 1));
                while k2 == k1
                    k2 = randi(size(pb, 1));
                end
                if pby(k1,1) < pby(k2,1)
                    add = k1;
                    update = k2;
                else
                    add = k2;
                    update = k1;
                end
                next_samplex = [next_samplex;pb(add,:)]; %add直接加入
                next_sampley = [next_sampley;pby(add,:)];
                next_v = [next_v;pbv(add,:)];
                pbv(update,:) = rand(1, num_vari) .* pbv(update,:) + rand(1, num_vari) .* (pb(add,:) - pb(update,:));    %update需要更新
                xl = pb(update,:) + pbv(update,:);
                xl(xl > upper_bound(func_num)) = upper_bound(func_num);
                xl(xl < lower_bound(func_num)) = lower_bound(func_num);% 范围检查
                if evalutions <= max_evalutions
                    yl = benchmark_func(xl,func_num);
                    evalutions = evalutions + 1;
                    if yl < fmin
                        fmin = yl;
                    end
                    arr_fmin = [arr_fmin;fmin];
                    fprintf("func_num:%d evalutions:%d fmin: %.4f\n", func_num, evalutions, fmin);
                    next_samplex = [next_samplex;xl];   %加到下一次迭代中
                    next_sampley = [next_sampley;yl];
                    next_v = [next_v;pbv(update,:)];
                else
                    break;
                end

                % 从pb中去掉k1和k2,从后面删除,防止后面删除的数索引发生变化
                if k1 > k2
                    first = k1;
                    second = k2;
                else
                    first = k2;
                    second = k1;
                end
                pb(first, :) = [];
                pby(first, :) = [];
                pbv(first, :) = [];
                pb(second, :) = [];
                pby(second, :) = [];
                pbv(second, :) = [];
            end
        end
        
        % 更新sample_x,v和sample_y
        if size(next_samplex,1) ~= num_initial   %代表程序结束
            break;
        else
            sample_x = next_samplex;
            sample_y = next_sampley;
            v = next_v;
        end
    end
    plot(arr_fmin);
end